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When a simple excitable system is continuously stimulated by a Poissonian external source, the 
response function (mean activity versus stimulus rate) generally shows a linear saturating shape. 
This is experimentally verified in some classes of sensory neurons, which accordingly present a small 
dynamic range (defined as the interval of stimulus intensity which can be appropriately coded by 
the mean activity of the excitable element), usually about one or two decades only. The brain, 
on the other hand, can handle a significantly broader range of stimulus intensity, and a collective 
phenomenon involving the interaction among excitable neurons has been suggested to account for 
the enhancement of the dynamic range. Since the role of the pattern of such interactions is still 
unclear, here we investigate the performance of a scale-free (SF) network topology in this dynamic 
range problem. Specifically, we study the transfer function of disordered SF networks of excitable 
Greenberg-Hastings cellular automata. We observe that the dynamic range is maximum when 
the coupling among the elements is critical, corroborating a general reasoning recently proposed. 
Although the maximum dynamic range yielded by general SF networks is slightly worse than that 
of random networks, for special SF networks which lack loops the enhancement of the dynamic 
range can be dramatic, reaching nearly five decades. In order to understand the role of loops on 
the transfer function we propose a simple model in which the density of loops in the network can 
be gradually increased, and show that this is accompanied by a gradual decrease of dynamic range. 

PACS numbers: 87.19.La, 87.18. Sn, 89. 75. He 89.20.-a 



I. INTRODUCTION 

Recent applications of tools from Statistical Physics 
have brought about new perspectives to theoretical Neu- 
roscience. On the one hand, networks of simplified neu- 
ron models seem to capture essential features of collective 
neuronal dynamics [J, [3, H, 0, Q , very often being also 
amenable to analytical calculations via mean field ap- 
proximations or Fokker-Planck equations [1,0, Si- On the 
other hand, experimental data from real neural networks 
have yielded extremely interesting results when analyzed 
within the framework of complex networks, often reveal- 
ing the small-world character for structural (i.e. anatom- 
ical) connectivity in different spatial scales 0, [H, |Tl[ , as 
well as scale-free characteristics for functional connectiv- 
ity both in in-vitro [T^ and in fMRI data (see 
for a recent review). 

In the context of modelling, one intriguing question 
in Neuroscience regards the stunning ability that brains 
have to cope with sensory stimuli that vary over many 
orders of magnitude Q • The experimental evidence sup- 
porting this claim has been accumulating for about a 
century in the Psychophysics literature: the perception 
of a given stimulus grows with a power law of the stim- 
ulus intensity (Stevens law), with an exponent (Stevens 
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exponent) which is typically < 1, implying low-stimulus 
amplification and large dynamic range [16[. This is in 
stark contrast with the poor performance of single neu- 
rons: as a function of the stimulus intensity, the mean 
firing rate of sensory neurons experimentally shows the 
linear saturating shape that one expects for general ex- 
citable systems, so their responses consistently have a 
small dynamic range, typically about one or two decades 
only [l3, [iH [l^, Hyl- How can these two experimen- 
tal results be reconciled? A solution which has been 
proposed for this apparent paradox involves a collective 
phenomenon. The idea is that if excitable elements with 
small dynamic range are coupled, signal propagation in 
the network amplifies the average activity, as compared 
to that of an isolated node. This collectively leads to a 
significant enhancement of dynamic range, thus provid- 
ing a possible solution to a problem faced by biological as 
well as artificial sensors: how to code for several orders 
of magnitude of stimulus intensity, starting from narrow- 
coding elements 0, i, i, i, 01 ■ 

The reasoning underlying the enhancement of dynamic 
range is very general and applies to essentially any net- 
work topology. Consider the limit of very weak stimulus, 
where each excitable element has a small probability of 
being excited. By coupling the elements, a single stimu- 
lus event will be amplified to neighboring sites, which will 
further amplify it, and so forth. If the coupling strength 
is small, this excitable wave will eventually die out, but 
the overall network activity (the response to the stim- 
ulus) will still be larger than that of the isolated ele- 
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ment that originally received the stimulus. The larger 
the coupling strength, the larger the amplification, and 
so the sensitivity and the dynamic range of the response 
curve initially increase with coupling. There is, however, 
a critical value of the coupling above which self-sustained 
activity becomes stable. Above this (typically second 
order) nonequilibrium phase transition, the response of 
the network for weak stimulus is hindered, because it 
is masked by the self-sustained activity of the network. 
This gets worse and worse as the coupling increases, so 
above criticality the dynamic range of the response curve 
decreases with increasing coupling strength. Therefore, 
the dynamic range is optimal at criticality [2]. 

The above mechanism has been tested in regular 0, 
H) 0, [1] as well as random Q networks of excitable el- 
ements. The maximum enhancement in dynamic range 
is about 100% in one-dimensional networks and 50% 
in random graphs It is not a priori clear how the 
performance depends on the network topology and, in 
particular, which one gives the best results. Given the 
potential applications of the mechanism to artificial sen- 
sors, as well as the relevance to Neuroscience, in this pa- 
per we study the performance of a scale-free topology in 
the dynamic range problem. We show that a particular 
class of scale-free networks, those with no loops, yield the 
best performance known so far. We investigate the role 
of loops on the dynamic range by introducing a slightly 
modified version of the Barabasi- Albert scale-free model 
where we can now tune the amount of loops in the net- 
work. 

The paper is organized as follows. In section |TT] we 
introduce the model and give a precise definition of the 
dynamic range. Results are discussed in section IIIII for 
the standard Barabasi- Albert scale-free model (jlll Ap as 
well as for a slightly modified "loop-diluted" version that 
we introduce pil Bl) . Our concluding remarks are pre- 
sented in section HVl 



II. THE MODEL 

We consider a variant of the Greenberg-Hastings cellu- 
lar automaton |2l|, which is one of the simplest models 
of an excitable system and can be used in large-scale sim- 
ulations. In the model, each excitable node i = 1, . . . , N 
could represent either a neuron, an active dendritic patch 
or even sub-cellular excitable processes. Each node can 
be in one of n states: Xi = is the quiescent state (e.g. 
a polarized neuron), Xi — 1 is the excited state (e.g. 
a spiking neuron) and Xi — 2, . . . , n — 1 are refractory 
states (e.g. a hyperpolarized neuron). Once a site is ex- 
cited {xi = 1), it deterministically goes through the next 
n — 2 refractory states, after which it jumps to the quies- 
cent state Xi = (the automaton is therefore cyclic [22|). 
Each node is independently excited by a stochastic exter- 
nal source, which mimics the effect of an stimulus on the 
lattice. We model the arrival of a suprathreshold stimu- 
lus by a Poisson process with rate r: at each time step r 



an attempt to stimulate a site occurs with probability 

A = 1 — exp(— rr) (1) 

(we adopt the arbitrary time scale of r = 1 ms, which 
is the characteristic time scale of a neuronal spike). We 
refer to the rate r as the stimulus intensity. In order to 
become excited in time t -f- r a given site has to be in 
state at time t. There are two different ways by which 
a site can be excited: by the continuous stimulation of 
the external source (with probability A per time step) 
or by stimulus propagation from its excited neighbors. 
Thus the probability that a quiescent site i is excited in 
the next time step is 

P,(i + r) = 1 - (1 - A) - V^i)K^Af)^ 1) : (2) 

where ^(a, 6) is the Kronecker delta, fc^ is the number of 
neighbors (connectivity or degree) of site i and is the 
probability that excitation from site j gets transmitted 
to site i. There is quenched disorder in the coupling: 
the probabilities are initially drawn from a uniform 
distribution in [0,2cr/i^] if 2a jK < 1, or [2<7/K - 1, 1] 
if 2a /K > 1, where K = (k) is the mean connectiv- 
ity of the network and a is the coupling parameter (for 
simplicity, we consider the case of bidirectional coupling 
Pij = Pji). Note that, in a mean field approximation, a 
coincides with the branching ratio, defined as the average 
of the number of descendant excitations divided by the 
number of ancestor excitations of each site. Such mean 
field approximation provides a quite satisfactory agree- 
ment with simulation results for random graph topologies 
(as expected) and accurately predicts a phase transition 
at = 1 1,0- 

The mean firing rate of the network is defined as 

F = T-^YI Pu where pt = iV"^ Ef 1) is the 
instantaneous density of active (excited) sites and T is 
a given time window for measurements (we have used 
N = 10**, T — 10"* steps and n = 5 states in most 
simulations). We refer to F{r) as the response function 
(or transfer function) of the network. It typically shows 
the sigmoidal shape in a log-linear scale exemplified in 
Fig. [Tfa), with a baseline activity Fq = lim^^o F{f) and 
saturation at F„iax = limr^oo F{r). The dynamic range 
A of the response function is defined as the width (mea- 
sured in dB) in stimulus intensity r which can be "ap- 
propriately coded" by F. In the biological literature, this 
is usually operationalized as follows [ItI . [Tsj : by letting 
Fx = Fq + x{Fmax - Fq) , where < a; < 1, and r^ be 
the corresponding stimulus intensity, {F(rx) = F^, see 
triangles in Fig.[lja) for an example), the dynamic range 
is defined as 

A=101ogio(^^) , (3) 
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therefore excluding stimuli whose response is just above 
baseline (r < ro.i) or too close to saturation (r > rg.g). 
For an isolated Greenberg-Hastings excitable node, one 
can easily show that the dynamic range is A < 19 dB 0, 



III. RESULTS 
A. Barabasi-Albert networks 

We consider scale-free networks [l^ of such excitable 
elements. Several investigations show that distinct 
systems such as World-Wide Web [23j], scientific [2^ . 
metapopulation dynamics [1^ and biochemical net- 
works [27I . [28l | self-organize into a scale- free configuration 
[2^, which means that the probability Pk that a given 
node has k edges follows a power-law distribution like 



Pi. cc k 



(4) 



Measurements in real systems estimate 7 in the range 
[2,3]. Eq. (HI basically means that poorly-connected 
nodes are most frequent in the network than well- 
connected nodes (hubs). 

To establish scale-free networks, we follow the stan- 
dard algorithm by Barabasi and Albert (B A) ^2§\ , which 
regards preferential attachment and growth as mecha- 
nisms for the emergence of the scale-free character. In 
this algorithm, the resulting networks display connectiv- 
ity distribution according to P^ cx k^^. The parameters 
of the BA model are the number of nodes N and to, 
which corresponds to the number of links that a newly 
introduced node adds to the network. These m links are 
most probably attached to those nodes with an already 
large number of edges. 

Figure [1] shows the results for to = 10. For small val- 
ues of cr, the response function F{r) increases linearly 
for weak stimulus. This linearity can be easily inter- 
preted: each stimulus arrival generates an excitable wave 
that will have a finite lifetime and will die before another 
wave is created. For stronger stimulus (larger r) linear- 
ity breaks down, since there is interaction among waves, 
which partially annihilate each other. For very large r, 
the firing rates reach a saturation value which scales with 
the inverse of the refractory period, F„iax — 
As the value of a increases, so does the lifetime of an ex- 
citable wave, leading to larger amplification of weak stim- 
uli and a corresponding enhancement of dynamic range 
(see Fig. [Hfor cr < 0.5). When cr = CTc, the lifetime of 
the excitable waves effectively diverges and the system 
undergoes a second order phase transition (notice the 
change in the exponent in the filled circles of Fig.[T](b)). 
For a > ac, any perturbation in the network will lead to 
a stable self sustained activity, Fo > (inset of Fig.[Hb)) 
which, as explained in section [H leads to smaller values 
of the dynamic range as the coupling increases [Tj (see 
Fig. m for cr > 0.5). 
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FIG. 1: (a) Response functions for BA scale-free networks 
with m — 10: mean firing rate F versus stimulus rate r. 
Different curves denote different values of the branching pa- 
rameter a: from bottom to top, a = 0.1, 0.3, . . . , 1.9. Filled 
circles: a — 0.5, which is close to the critical value. The hori- 
zontal lines exemplify how the dynamic range A is calculated 
for (T = 1.9 (filled triangles), (b) Log- log version of (a). The 
dashed line shows an exponent 1/2. Inset: Self-sustained ac- 
tivity Fo versus cr, illustrating the phase transition close to 

(Tc ~ 0.5. 



One observes that, differently from random graph 
topologies 0, the transition for scale-free excitable net- 
works occurs at crc < 1. We speculate that this is due to 
the hubs, which have a local branching ratio cr^ — p^j 
larger than unit even for ct < 1 and could therefore fa- 
cilitate the transition. It is also interesting to note that 
deviations from mean field behavior have been predicted 
for the contact process (CP) in a scale-free network [30l |. 
Apart from the refractory period and the disorder, the 
CP is similar to the model we study here (in the sense 
that it has a unique absorbing state with no symmetries). 
In Fig. [1] however, the response exponent at criticality 
(defined by F{r;ac) ^ r^^^'^) seems to be compatible 
with the mean field value l/Sh = 1/2 [22]. 

Results in Fig. [T] are typical, similar curves are ob- 
tained for any m > 1. The performance of these scale- 
free networks in enhancing the dynamic range is poor: 
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FIG. 2: Dynamic range A versus branching parameter a for 
distinct network topologies: BA scale-free networks (open cir- 
cles) and Erdos-Renyi random graphs (filled circles) with ap- 
proximately the same mean connectivity. Inset: Response 
function for BA scale-free networks with m = 1 and a — 2 for 
N = 10* (filled circles) and = 3 x 10* (open circles). The 
solid line shows a slope ~ 0.07. 



while the dynamic range of isolated excitable elements is 
A(a = 0) ~ 16.7 dB, the network (optimal) performance 
at criticality is only A(crc) — 20.8 dB, an enhancement 
of less than a decade. This is slightly worse than the 
enhancement produced by random networks with equiv- 
alent size and average connectivity, as can be seen in the 
curves A{a) of Fig. [2l 

The case m = 1, however, is particularly interesting. 
Notice that in this situation the network is still scale- free, 
but does not comprise any loop in its structure and conse- 
quently has a tree-like pattern. This condition, together 
with the deterministic nature of each excitable node after 
excitation, prevents the phase transition to self-sustained 
activity from occurring [Ij], a fact that has also been ob- 
served in one-dimensional excitable networks @ . In these 
conditions the only transition occurs at cr = amax = K 
(= 2m for scale-free networks), whereby propagation of 
excitable waves becomes deterministic (ballistic) . There- 
fore low-stimulus amplification increases steadily with cr, 
but in the absence of self-sustained activity (Fg = 0). 
This allows the dynamic range to increase monotonically 
with cr, reaching values near 50 dB, which is the largest 
value obtained so far in excitable network models. 



B. Loop-diluted model 

As we observe a remarkable difference between the dy- 
namic range of scale-free networks with m — 1 and other 
values of to, and the former has a typical feature (non- 
existence of loops) which is not present in to > 1 topolo- 
gies, we are interested in investigating the role of loops in 
the response functions of the networks. For this purpose, 
we propose a variant of the BA model which is referred 
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FIG. 3: Occurrences of 3- and 4-site patterns in the loop- 
diluted scale-free networks (filled symbols) and their random- 
ized counterparts (open symbols), (a) Number of triangles 
versus p and (b) Number of squares versus p. Points (bars) 
represent the mean (standard deviation) over 5 realizations 
for A'^ = 1000. Patterns are sequentially numbered (1-8) for 
further reference in the text. 



to as loop-diluted model. In the model each new node 
is added according to the usual preferential attachment 
rule, but can have to = 1 or to = 2 links according to the 
probability distribution 



P(to) = (1 - p)5,n,i + P6m,2 



(5) 



where p is the probability of having two edges. So p 
adjusts the amount of loops in the network and the case 
p = recovers the structure with no loops. The mean 
degree is now K — 2 (to) = 2(1 + p). 

Notice that, since all sites belong to a single giant 
component, the average number of loops created by each 
newly added site is bounded from below by p. Each new 
two-edged node can give rise to loops of any size, but 
the relationship between parameter p and the number of 
loops becomes already apparent in a simple 3-site motif 
analysis 31J . Figure [3ja) shows the mean number of tri- 
angles as a function of p (calculated by the free software 
available at www.weizman.ac.il/mcb/UriAlon). As ex- 
pected, this is a monotonically increasing function, which 



5 




FIG. 4: Dynamic range A versus branching parameter cr for 
the loop-diluted model and different values of the probability 
p that a newly added node has two incoming edges. Inset; 
maximum value of the dynamic range as a function of p. 

nevertheless stays well below Np, hinting that most loops 
comprise more than three sites. The same qualitative 
scenario is observed when we plot the number of squares 
(Fig-Etb)), which are more abundant than triangles. In 
both cases, we notice that for equally sized randomized 
graphs (which preserve the degrees of every node [sH). 
the numbers of triangles and squares are considerably 
larger. This means that triangles and squares are actu- 
ally anti-motifs in the loop-diluted model [s^ ]. In fact, 
this is true for all patterns that contain loops (numbered 
2, 5, 6 and 8 in Fig. [3]), except for pattern 7, which cannot 
occur according to the growth rules of the model. The 
occurrences of patterns 1, 3 and 4 in the loop-diluted 
model and in the randomized networks are statistically 
indistinguishable (therefore they are neither motifs nor 
anti- motifs - data not shown). 

Figure |4] displays the dynamic range A as a function 
of the coupling a for some values of p. From the figure, 
we clearly notice that the insertion of loops by increasing 
the probability p has a striking effect on the dynamic 
range. This effect is already mensurable for values of 
p such that pN ^ 1. The peak value of the dynamic 
range Amaxip) = maxo-A(a-;p) = A{ac{p);p) seems to 
decrease logarithmically with p. 

IV. CONCLUDING REMARKS 

Recently, some investigations have addressed the en- 
largement in average activity of excitable elements by 
coupling these entities and so giving form to a new larger 
and more sensitive unit. Although this collective phe- 
nomenon has been widely accepted, little is known about 



the way the arrangement of connections among the ex- 
citable elements acts physically on the system dynamics. 
The creation of more robust and functional units from 
smaller units (whose pattern of interactions is a deter- 
mining aspect) is of course not exclusive of Neuroscience. 
For instance, the hypercycle, a catalytic feedback net- 
work whereby each element helps the replication of the 
next one in a regulatory cycle closing on itself, has been 
pondered as an alternative resolution for the information 
crisis in prebiotic evolution [sl, [13] . We believe that all 
the recent contributions on this issue have a bearing on 
a more general context, that is, the understanding of the 
interplay between system dynamics and the underlying 
interaction network topologies. We hope that our con- 
tribution gives a small step in this direction, when we 
corroborate that the amount of loops in network struc- 
ture could be a key topological feature. 

We have presented simulation results for the transfer 
function of excitable scale-free networks. The behavior 
of the dynamic range A as a function of the coupling 
a shows the general behavior predicted in Ref. [7|: in 
the subcritical regime (cr < CTc) A.{a) increases, while in 
the supercritical regime {a > cjc) A((t) decreases. The 
maximum value is obtained at criticality, but for scale- 
free networks with m > 1 this result is even smaller than 
that for a random graph. 

For m = 1 the phase transition to self-sustained activ- 
ity disappears, and the dynamic range increases steadily, 
reaching its maximum value when excitable waves be- 
come deterministic. This suggests that the presence of 
loops in the network could be a relevant feature in deter- 
mining the dynamic range of its transfer function. We 
have introduced a simple extension to the BA scale-free 
model which allows one to interpolate between m — 1 
and m = 2, showing that dynamic range decreases as the 
density of loops increases. This reinforces the need to 
study other topologies with tree structure, which abound 
in biological structures. 

It remains at present unclear whether the response ex- 
ponent for m > 1 is indeed compatible with the mean 
field universality class (even though this seems to be sup- 
ported by recent simulation results in Ref. (ssj , which 
independently addressed a similar problem). Also, for 
TO = 1 at maximum coupling, the response function 
seems to be governed by a power law with a much smaller 
exponent (see inset of Fig. [2]) which might not belong 
to the directed percolation universality class. We be- 
lieve that a detailed study of the critical exponents of 
excitable scale-free networks is still lacking and should 
be dealt with in the future. 
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